/** Author: David Powell
This file produces Figures E.1-E.3
**/


clear all
set more off
set mat 800
set seed 8721

global dir "/jules/b/dpowell"
global DATA "${dir}/purdue/replication/DATA"
global OUTPUT "${dir}/purdue/replication/output"


use ${DATA}/finaldata if year>1990

gen overdose_rate=100000*(overdose/totpop)
gen opioid_rate=100000*opioids/totpop


gen post1=(year>1995&year<2001)*nontripl
gen post2=(year>2000&year<2011)*nontripl
gen post3=(year>2010)*nontripl


tab stfips, gen(ss) nof
tab year, gen(ttt) nof
drop ttt13
foreach nnn of numlist 1983/2017 {
	gen tt_`nnn'=nontripl*(year==`nnn')
}


**saving results into these matrices
matrix BBB=J(10000,2,.)
matrix CCC=J(10000,6,.)
matrix DDD=J(10000,2,.)
matrix TTT=J(10000,6,.)


set seed 30985231



***get values when treatment is correctly assigned***
reg overdose_rate ss* tt_1995 ttt*  [aw=totpop] if year==1991 | year==1995, cluster(stf)
matrix BBB[1,1]=_b[tt_1995]
matrix BBB[1,2]=_b[tt_1995]/_se[tt_1995]

reg opioid_rate ss* tt_1995 ttt*  [aw=totpop] if year==1991 | year==1995, cluster(stf)
matrix DDD[1,1]=_b[tt_1995]
matrix DDD[1,2]=_b[tt_1995]/_se[tt_1995]


reg overdose_rate ss* post1 post2 post3 ttt*  [aw=totpop] if year>=1991, cluster(stf)
matrix CCC[1,1]=_b[post1]
matrix CCC[1,4]=_b[post1]/_se[post1]

matrix CCC[1,2]=_b[post2]
matrix CCC[1,5]=_b[post2]/_se[post2]

matrix CCC[1,3]=_b[post3]
matrix CCC[1,6]=_b[post3]/_se[post3]



reg opioid_rate ss* post1 post2 post3 ttt*  [aw=totpop] if year>=1991, cluster(stf)
matrix TTT[1,1]=_b[post1]
matrix TTT[1,4]=_b[post1]/_se[post1]

matrix TTT[1,2]=_b[post2]
matrix TTT[1,5]=_b[post2]/_se[post2]

matrix TTT[1,3]=_b[post3]
matrix TTT[1,6]=_b[post3]/_se[post3]



***get placebo values***
local i=2
gen oldnontripl=nontripl
qui while `i'<=10000 {
	drop nontripl  post*

	**random assignment of non-triplicate status	
	gen tmp=uniform() if year==1991 
	egen tmp2=rank(tmp), unique by(year oldnontripl)
	egen tmp3=max(tmp2), by(stf)
	
	gen nontripl=(tmp3<=41 | oldnontripl==0)

	gen post1=nontripl*(year>=1996 & year<=2000)
	gen post2=nontripl*(year>=2001 & year<=2010)
	gen post3=nontripl*(year>=2011)
	replace tt_1995=nontripl*(year==1995)

	drop tmp*

	reg overdose_rate ss* post1-post3 ttt*  [aw=totpop] if year>=1991, cluster(stf)
	matrix CCC[`i',1]=_b[post1]
	matrix CCC[`i',4]=_b[post1]/_se[post1]
	matrix CCC[`i',2]=_b[post2]
	matrix CCC[`i',5]=_b[post2]/_se[post2]
	matrix CCC[`i',3]=_b[post3]
	matrix CCC[`i',6]=_b[post3]/_se[post3]

	reg opioid_rate ss* post1-post3  ttt*  [aw=totpop] if year>=1991, cluster(stf)
	matrix TTT[`i',1]=_b[post1]
	matrix TTT[`i',4]=_b[post1]/_se[post1]
	matrix TTT[`i',2]=_b[post2]
	matrix TTT[`i',5]=_b[post2]/_se[post2]
	matrix TTT[`i',3]=_b[post3]
	matrix TTT[`i',6]=_b[post3]/_se[post3]



	reg overdose_rate ss* tt_1995 ttt*  [aw=totpop] if year==1991 | year==1995, cluster(stf)
	matrix BBB[`i',1]=_b[tt_1995]
	matrix BBB[`i',2]=_b[tt_1995]/_se[tt_1995]

	reg opioid_rate ss* tt_1995 ttt*  [aw=totpop] if year==1991 | year==1995, cluster(stf)
	matrix DDD[`i',1]=_b[tt_1995]
	matrix DDD[`i',2]=_b[tt_1995]/_se[tt_1995]


	local i=`i'+1
}

svmat CCC
svmat TTT
svmat BBB
svmat DDD


sum BBB1 if _n==1
local rrr=r(mean)
centile BBB1, centile(97.5)
local rrr1=r(c_1)
centile BBB1, centile(2.5)
local rrr0=r(c_1)

histogram BBB1 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick) lstyle(foreground)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) xtitle("") xlabel(-5(1)2) fraction 
graph export figE3A.eps, replace
graph export figE3A.pdf, replace

# delimit cr

sum BBB2 if _n==1
local rrr=r(mean)
centile BBB2, centile(97.5)
local rrr1=r(c_1)
centile BBB2, centile(2.5)
local rrr0=r(c_1)

histogram BBB2 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)  lstyle(foreground)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick))  xtitle("") xlabel(-6(1)5) fraction 
graph export figE3B.eps, replace
graph export figE3B.pdf, replace


sum DDD1 if _n==1
local rrr=r(mean)
centile DDD1, centile(97.5)
local rrr1=r(c_1)
centile DDD1, centile(2.5)
local rrr0=r(c_1)

histogram DDD1 if _n>1 & _n<=10000,  graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)  lstyle(foreground)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick))  xtitle("") xlabel(-2(1)1) fraction 
graph export figE3C.eps, replace
graph export figE3C.pdf, replace


sum DDD2 if _n==1
local rrr=r(mean)
centile DDD2, centile(97.5)
local rrr1=r(c_1)
centile DDD2, centile(2.5)
local rrr0=r(c_1)
sum DDD2

histogram DDD2 if _n>1 & _n<=10000 & DDD2>-6, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)  lstyle(foreground)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) xtitle("") xlabel(-6(1)6) fraction 
graph export figE3D.eps, replace
graph export figE3D.pdf, replace



sum CCC1 if _n==1
local rrr=r(mean)
centile CCC1, centile(97.5)
local rrr1=r(c_1)
centile CCC1, centile(2.5)
local rrr0=r(c_1)

histogram CCC1 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist1_noxx.gph, replace) xtitle("") xlabel(-2(.5)1.5) fraction 
graph export figE1A.eps, replace
graph export figE1A.pdf, replace


sum CCC2 if _n==1
local rrr=r(mean)
centile CCC2, centile(97.5)
local rrr1=r(c_1)
centile CCC2, centile(2.5)
local rrr0=r(c_1)

histogram CCC2 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist2_noxx.gph, replace) xtitle("") fraction  xlabel(-5(1)5)
graph export figE1B.eps, replace
graph export figE1B.pdf, replace


sum CCC3 if _n==1
local rrr=r(mean)
centile CCC3, centile(97.5)
local rrr1=r(c_1)
centile CCC3, centile(2.5)
local rrr0=r(c_1)
histogram CCC3 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist3_noxx.gph, replace) xtitle("") fraction xlabel(-12(2)8)
graph export figE1C.eps, replace
graph export figE1C.pdf, replace


sum CCC4 if _n==1
local rrr=r(mean)
centile CCC4, centile(97.5)
local rrr1=r(c_1)
centile CCC4, centile(2.5)
local rrr0=r(c_1)

histogram CCC4 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist4_noxx.gph, replace) xtitle("") fraction xlabel(-6(1)4)
graph export figE2A.eps, replace
graph export figE2A.pdf, replace

sum CCC5 if _n==1
local rrr=r(mean)
centile CCC5, centile(97.5)
local rrr1=r(c_1)
centile CCC5, centile(2.5)
local rrr0=r(c_1)

histogram CCC5 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist5_noxx.gph, replace) xtitle("") fraction xlabel(-6(1)5.5)
graph export figE2B.eps, replace
graph export figE2B.pdf, replace

sum CCC6
sum CCC6 if _n==1
local rrr=r(mean)
centile CCC6, centile(97.5)
local rrr1=r(c_1)
centile CCC6, centile(2.5)
local rrr0=r(c_1)
histogram CCC6 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) xlabel(-6(1)6) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist6_noxx.gph, replace) xtitle("") fraction xlabel(-10(1)6)
graph export figE2C.eps, replace
graph export figE2C.pdf, replace


sum TTT1 if _n==1
local rrr=r(mean)
centile TTT1, centile(97.5)
local rrr1=r(c_1)
centile TTT1, centile(2.5)
local rrr0=r(c_1)

histogram TTT1 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist1a_noxx.gph, replace) xtitle("") xlabel(-2.5(.5)1) fraction 
graph export figE1D.eps, replace
graph export figE1D.pdf, replace

sum TTT2 if _n==1
local rrr=r(mean)
centile TTT2, centile(97.5)
local rrr1=r(c_1)
centile TTT2, centile(2.5)
local rrr0=r(c_1)

histogram TTT2 if _n>1 & _n<=10000,  graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist2a_noxx.gph, replace) xtitle("") fraction  xlabel(-6(1)3)
graph export figE1E.eps, replace
graph export figE1E.pdf, replace

sum TTT3 if _n==1
local rrr=r(mean)
centile TTT3, centile(97.5)
local rrr1=r(c_1)
centile TTT3, centile(2.5)
local rrr0=r(c_1)
histogram TTT3 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist3a_noxx.gph, replace) xtitle("") fraction xlabel(-12(2)6)
graph export figE1F.eps, replace
graph export figE1F.pdf, replace


sum TTT4 if _n==1
local rrr=r(mean)
centile TTT4, centile(97.5)
local rrr1=r(c_1)
centile TTT4, centile(2.5)
local rrr0=r(c_1)

histogram TTT4 if _n>1 & _n<=10000,  graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist4a_noxx.gph, replace) xtitle("") fraction xlabel(-7(1)4)
graph export figE2D.eps, replace
graph export figE2D.pdf, replace

sum TTT5 if _n==1
local rrr=r(mean)
centile TTT5, centile(97.5)
local rrr1=r(c_1)
centile TTT5, centile(2.5)
local rrr0=r(c_1)

histogram TTT5 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist5a_noxx.gph, replace) xtitle("") fraction xlabel(-6(1)5)
graph export figE2E.eps, replace
graph export figE2E.pdf, replace


sum TTT6
sum TTT6 if _n==1
local rrr=r(mean)
centile TTT6, centile(97.5)
local rrr1=r(c_1)
centile TTT6, centile(2.5)
local rrr0=r(c_1)
histogram TTT6 if _n>1 & _n<=10000, graphregion(color(white)) xline(`rrr', lcolor(blue) lwidth(thick)) ylabel(,nogrid) xline(`rrr0', lpattern(dash) lwidth(thick)) xline(`rrr1', lpattern(dash) lwidth(thick)) saving(hist6a_noxx.gph, replace) xtitle("") fraction xlabel(-9(1)4)
graph export figE2F.eps, replace
graph export figE2F.pdf, replace





foreach var of varlist CCC* TTT* BBB* DDD* {
	egen rk_`var'=rank(`var'), field
	egen rk1_`var'=rank(abs(`var')), field
}

***MAIN RESULTS***
sum rk* if _n==1

for num 1/6: gen sAAAX=CCCX[1]
for num 1/2: gen sBBBX=BBBX[1]
for num 1/2: gen sDDDX=DDDX[1]
for num 1/6: gen SSSX=TTTX[1]




***joint tests
count if sAAA1 < CCC1 & sAAA2 < CCC2 & sAAA3 < CCC3 & _n>1 & CCC1<.
count if abs(sAAA1) < abs(CCC1) & abs(sAAA2) < abs(CCC2) & abs(sAAA3) < abs(CCC3) & _n>1 & CCC1<.

count if sAAA4 < CCC4 & sAAA5 < CCC5 & sAAA6 < CCC6 & _n>1 & CCC1<.
count if abs(sAAA4) < abs(CCC4) & abs(sAAA5) < abs(CCC5) & abs(sAAA6) < abs(CCC6) & _n>1 & CCC1<.

count if SSS1 < TTT1 & SSS2 < TTT2 & SSS3 < TTT3 & _n>1   & CCC1<.
count if abs(SSS1) < abs(TTT1) & abs(SSS2) < abs(TTT2) & abs(SSS3) < abs(TTT3) & _n>1   & CCC1<.

count if SSS4 < TTT4 & SSS5 < TTT5 & SSS6 < TTT6 & _n>1   & CCC1<.
count if abs(SSS4) < abs(TTT4) & abs(SSS5) < abs(TTT5) & abs(SSS6) < abs(TTT6) & _n>1   & CCC1<.
